This is a documentation for analyses of scATAC-seq data, generated from rat metrial gland tissues on gestational day (GD) 15.5 and 19.5.
library(Seurat)
library(Signac)
library(ggplot2)
library(GenomicRanges)
set.seed(1234)
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/rnor6.rda")
colnames(rnor6)[1] <- "gene_id"
rnor6.ranges <- makeGRangesFromDataFrame(rnor6, seqnames.field = "chromosome_name",
start.field = "start_position",
end.field = "end_position",
strand.field = "strand",
keep.extra.columns = TRUE)
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/scRNA-seqFiles/gd15.5_res0.8.rda")
data <- RenameIdents(data, `0` = "Other cells", `1` = "Other cells",
`2` = "Other cells", `3` = "Other cells", `4` = "Other cells",
`5` = "Other cells", `6` = "Natural killer cells",
`7` = "Macrophage cells", `8` = "Macrophage cells",
`9` = "Natural killer cells", `10` = "Endothelial cells",
`11` = "Macrophage cells", `12` = "Other cells",
`13` = "Endothelial cells", `14` = "Smooth muscle cells",
`15` = "Endothelial cells", `16` = "Other cells",
`17` = "Smooth muscle cells", `18` = "Other cells",
`19` = "Natural killer cells", `20` = "Other cells",
`21` = "Other cells", `22` = "Other cells",
`23` = "Invasive trophoblast cells", `24` = "Other cells",
`25` = "Other cells", `26` = "Smooth muscle cells",
`27` = "Other cells", `28` = "Other cells")
Idents(data) <- factor(Idents(data), levels = c("Other cells", "Macrophage cells", "Smooth muscle cells",
"Endothelial cells", "Natural killer cells", "Invasive trophoblast cells"))
rats.RNA <- data
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/2_mergedSamples/2_MACS2Peaks/GD15.5/15.5_1-2-3-merged.rda")
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/rnor6.ranges.rda")
rats <- unintegrated
DefaultAssay(rats) <- "combined"
mcols(rnor6.ranges, level="within")[, "gene_id"] <- mcols(rnor6.ranges, level="within")[, "ensembl_gene_id"]
Annotation(rats) <- rnor6.ranges
gene.activities <- GeneActivity(rats)
rats[['RNA']] <- CreateAssayObject(counts = gene.activities[, colnames(rats)])
#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')
rats <- NormalizeData(object = rats, assay = 'RNA', normalization.method = 'LogNormalize', scale.factor = median(rats$nCount_RNA))
DefaultAssay(rats) <- 'RNA'
transfer.anchors <- FindTransferAnchors(reference = rats.RNA, query = rats, reduction = 'cca')
predicted.labels <- TransferData(anchorset = transfer.anchors, refdata = Idents(rats.RNA), weight.reduction = rats[['lsi']], dims=2:20)
rats <- AddMetaData(object = rats, metadata= predicted.labels)
plot1 <- DimPlot(object = rats.RNA, cols = c("grey70", "orange3", "#718748", "#183D61", "#79704C", "red3")) + ggtitle('GD15.5 scRNA-seq')
plot2 <- DimPlot(object=rats, group.by = 'predicted.id', label = T, repel = T) + ggtitle('GD15.5 scATAC-seq')
plot1 + plot2
#save(rats, file="/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/gd15.5-ATACwithRNAlables.rda")
table(rats$predicted.id)
#>
#> Endothelial cells Invasive trophoblast cells
#> 3250 128
#> Macrophage cells Natural killer cells
#> 2682 1443
#> Other cells Smooth muscle cells
#> 16700 1118
Correlation check:
library(Seurat)
library(Signac)
library(Seurat)
library(S4Vectors)
library(GenomicRanges)
library(patchwork)
library(biomaRt)
library(ggplot2)
library(ggpubr)
#> Warning: package 'ggpubr' was built under R version 4.0.3
library(BSgenome.Rnorvegicus.Ensembl.rn6)
set.seed(1234)
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd15.5-ATACwithRNAlables.rda")
Idents(rats) <- rats@meta.data$predicted.id # change from old identities to those predicted by scRNA-seq
DefaultAssay(rats) <- 'MACSpeaks'
avgATAC <- AverageExpression(rats, assays = "RNA")
avgATAC <- as.data.frame(avgATAC)
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/scRNA-seqFiles/gd15.5_res0.8.rda")
avgRNA <- AverageExpression(data, assays = "RNA")
avgRNA <- as.data.frame(avgRNA)
Trophoblast:
avgRNA1 <- avgRNA$RNA.23
names(avgRNA1) <- rownames(avgRNA)
avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Invasive.trophoblast.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]
cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#>
#> Spearman's rank correlation rho
#>
#> data: avgRNA1 and avgATAC1
#> S = 4.4498e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.4378357
df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-15.5-tbcCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Invasive Trophoblast Cells at gd 15.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")
#dev.off()
Endothelial:
avgRNA1 <- rowMeans(avgRNA[, c("RNA.10", "RNA.13", "RNA.15")])
names(avgRNA1) <- rownames(avgRNA)
avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Endothelial.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]
cor.test(log10(avgRNA1+10^(-5)), log10(avgATAC1+10^(-5)), method = "spearman", exact = FALSE)
#>
#> Spearman's rank correlation rho
#>
#> data: log10(avgRNA1 + 10^(-5)) and log10(avgATAC1 + 10^(-5))
#> S = 3.9856e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.4964745
df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-15.5-endoCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=FALSE) + ggtitle("Correlation in Endothelial Cells at gd 15.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")
#dev.off()
Macrophage:
avgRNA1 <- rowMeans(avgRNA[, c("RNA.7", "RNA.8", "RNA.11")])
names(avgRNA1) <- rownames(avgRNA)
avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Macrophage.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]
cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#>
#> Spearman's rank correlation rho
#>
#> data: avgRNA1 and avgATAC1
#> S = 3.6982e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.5327817
df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-15.5-macroCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Macrophage Cells at gd 15.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")
#dev.off()
NK:
avgRNA1 <- rowMeans(avgRNA[, c("RNA.6", "RNA.9", "RNA.19")])
names(avgRNA1) <- rownames(avgRNA)
avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Natural.killer.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]
cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#>
#> Spearman's rank correlation rho
#>
#> data: avgRNA1 and avgATAC1
#> S = 4.0055e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.4939641
df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-15.5-nkCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Natural Killer Cells at gd 15.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")
#dev.off()
Smooth muscle:
avgRNA1 <- rowMeans(avgRNA[, c("RNA.14", "RNA.17", "RNA.26")])
names(avgRNA1) <- rownames(avgRNA)
avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Smooth.muscle.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]
cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#>
#> Spearman's rank correlation rho
#>
#> data: avgRNA1 and avgATAC1
#> S = 4.382e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.4464026
df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-15.5-smCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Smooth Muscle Cells at gd 15.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")
#dev.off()
set.seed(1234)
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/scRNA-seqFiles/gd19.5-4-5-6-7_res0.8.rda")
data <- data2
data <- RenameIdents(data, `0` = "Other cells", `1` = "Macrophage cells",
`2` = "Other cells", `3` = "Other cells", `4` = "Other cells",
`5` = "Other cells", `6` = "Invasive trophoblast cells",
`7` = "Macrophage cells", `8` = "Smooth muscle cells",
`9` = "Natural killer cells", `10` = "Endothelial cells",
`11` = "Macrophage cells", `12` = "Other cells",
`13` = "Other cells", `14` = "Macrophage cells",
`15` = "Other cells", `16` = "Other cells",
`17` = "Macrophage cells", `18` = "Macrophage cells",
`19` = "Other cells", `20` = "Endothelial cells",
`21` = "Other cells", `22` = "Other cells",
`23` = "Other cells", `24` = "Other cells",
`25` = "Other cells", `26` = "Smooth muscle cells")
Idents(data) <- factor(Idents(data), levels = c("Other cells", "Macrophage cells", "Smooth muscle cells",
"Endothelial cells", "Natural killer cells", "Invasive trophoblast cells"))
rats.RNA <- data
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/2_mergedSamples/2_MACS2Peaks/GD19.5/19.5_4-5-6-merged.rda")
rats <- unintegrated
DefaultAssay(rats) <- "combined"
mcols(rnor6.ranges, level="within")[, "gene_id"] <- mcols(rnor6.ranges, level="within")[, "ensembl_gene_id"]
Annotation(rats) <- rnor6.ranges
gene.activities <- GeneActivity(rats)
rats[['RNA']] <- CreateAssayObject(counts = gene.activities[, colnames(rats)])
#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')
rats <- NormalizeData(object = rats, assay = 'RNA', normalization.method = 'LogNormalize', scale.factor = median(rats$nCount_RNA))
DefaultAssay(rats) <- 'RNA'
transfer.anchors <- FindTransferAnchors(reference = rats.RNA, query = rats, reduction = 'cca')
#> Warning in RunCCA.Seurat(object1 = reference, object2 = query, features =
#> features, : Running CCA on different assays
predicted.labels <- TransferData(anchorset = transfer.anchors, refdata = Idents(rats.RNA), weight.reduction = rats[['lsi']], dims=2:10)
rats <- AddMetaData(object = rats, metadata = predicted.labels)
plot1 <- DimPlot(object = rats.RNA, cols = c("grey70", "orange3", "#718748", "#183D61", "#79704C", "red3")) + ggtitle('GD19.5 scRNA-seq')
plot2 <- DimPlot(object=rats, group.by = 'predicted.id', label = T, repel = T) + ggtitle('GD19.5 scATAC-seq')
plot1 + plot2
#save(rats, file="/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/gd19.5-ATACwithRNAlables.rda")
table(rats$predicted.id)
#>
#> Endothelial cells Invasive trophoblast cells
#> 894 684
#> Macrophage cells Natural killer cells
#> 1265 223
#> Other cells Smooth muscle cells
#> 11125 197
Correlation check:
set.seed(1234)
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-ATACwithRNAlables.rda")
Idents(rats) <- rats@meta.data$predicted.id # change from old identities to those predicted by scRNA-seq
DefaultAssay(rats) <- 'MACSpeaks'
avgATAC <- AverageExpression(rats, assays = "RNA")
avgATAC <- as.data.frame(avgATAC)
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/scRNA-seqFiles/gd19.5-4-5-6-7_res0.8.rda")
avgRNA <- AverageExpression(data2, assays = "RNA")
avgRNA <- as.data.frame(avgRNA)
Trophoblast:
avgRNA1 <- avgRNA$RNA.6
names(avgRNA1) <- rownames(avgRNA)
avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Invasive.trophoblast.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]
cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#>
#> Spearman's rank correlation rho
#>
#> data: avgRNA1 and avgATAC1
#> S = 4.109e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.5258185
df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-19.5-tbcCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho") + ggtitle("Correlation in Invasive Trophoblast Cells at gd 19.5")
#dev.off()
Endothelial:
avgRNA1 <- rowMeans(avgRNA[, c("RNA.10", "RNA.20")])
names(avgRNA1) <- rownames(avgRNA)
avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Endothelial.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]
cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#>
#> Spearman's rank correlation rho
#>
#> data: avgRNA1 and avgATAC1
#> S = 3.9959e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.5388751
df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-19.5-endoCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Endothelial Cells at gd 19.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")
#dev.off()
Macrophage:
avgRNA1 <- rowMeans(avgRNA[, c("RNA.1", "RNA.7", "RNA.11", "RNA.17", "RNA.18")])
names(avgRNA1) <- rownames(avgRNA)
avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Macrophage.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]
cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#>
#> Spearman's rank correlation rho
#>
#> data: avgRNA1 and avgATAC1
#> S = 4.0517e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.532435
df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-19.5-macroCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Macrophage Cells at gd 19.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")
#dev.off()
NK:
avgRNA1 <- avgRNA[, c("RNA.9")]
names(avgRNA1) <- rownames(avgRNA)
avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Natural.killer.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]
cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#>
#> Spearman's rank correlation rho
#>
#> data: avgRNA1 and avgATAC1
#> S = 4.0328e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.5346098
df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-19.5-nkCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Natural Killer Cells at gd 19.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")
#dev.off()
Smooth muscle:
avgRNA1 <- rowMeans(avgRNA[, c("RNA.8", "RNA.26")])
names(avgRNA1) <- rownames(avgRNA)
avgATAC1 <- avgATAC[intersect(rownames(avgATAC), names(avgRNA1)), "RNA.Smooth.muscle.cells"]
names(avgATAC1) <- intersect(rownames(avgATAC), names(avgRNA1))
avgATAC1 <- avgATAC1[order(names(avgATAC1))]
avgRNA1 <- avgRNA1[intersect(rownames(avgATAC), names(avgRNA1))]
avgRNA1 <- avgRNA1[order(names(avgRNA1))]
cor.test(avgRNA1, avgATAC1, method = "spearman", exact = F)
#>
#> Spearman's rank correlation rho
#>
#> data: avgRNA1 and avgATAC1
#> S = 4.438e+11, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.4878489
df <- data.frame(avgRNA1=avgRNA1, avgATAC1=avgATAC1)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS4-19.5-smCorr.#pdf", height = 5, width = 6)
ggplot(df, aes(x=log10(avgRNA1+10^(-5)), y=log10(avgATAC1+10^(-5)))) + geom_point(alpha = 1/10) + geom_smooth(method=lm, se=FALSE, fullrange=TRUE) + ggtitle("Correlation in Smooth Muscle Cells at gd 19.5") + stat_cor(method = "spearman", label.y = -4, cor.coef.name = "rho")
#dev.off()
sessionInfo()
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> attached base packages:
#> [1] parallel stats4 stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] BSgenome.Rnorvegicus.Ensembl.rn6_0.0.1
#> [2] BSgenome_1.56.0
#> [3] rtracklayer_1.48.0
#> [4] Biostrings_2.56.0
#> [5] XVector_0.28.0
#> [6] ggpubr_0.4.0
#> [7] biomaRt_2.44.4
#> [8] patchwork_1.0.1
#> [9] GenomicRanges_1.40.0
#> [10] GenomeInfoDb_1.24.2
#> [11] IRanges_2.22.2
#> [12] S4Vectors_0.26.1
#> [13] BiocGenerics_0.34.0
#> [14] ggplot2_3.3.6
#> [15] Signac_1.5.0
#> [16] SeuratObject_4.0.4
#> [17] Seurat_4.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] utf8_1.2.2 reticulate_1.16 tidyselect_1.1.2
#> [4] RSQLite_2.2.1 AnnotationDbi_1.50.3 htmlwidgets_1.5.2
#> [7] grid_4.0.2 docopt_0.7.1 BiocParallel_1.22.0
#> [10] Rtsne_0.15 munsell_0.5.0 codetools_0.2-16
#> [13] ica_1.0-2 future_1.19.1 miniUI_0.1.1.1
#> [16] withr_2.5.0 colorspace_2.0-3 Biobase_2.48.0
#> [19] highr_0.9 knitr_1.39 rstudioapi_0.13
#> [22] ROCR_1.0-11 ggsignif_0.6.0 tensor_1.5
#> [25] listenv_0.8.0 labeling_0.4.2 slam_0.1-47
#> [28] GenomeInfoDbData_1.2.3 polyclip_1.10-0 bit64_4.0.5
#> [31] farver_2.1.0 vctrs_0.4.1 generics_0.1.2
#> [34] xfun_0.31 BiocFileCache_1.12.1 lsa_0.73.2
#> [37] ggseqlogo_0.1 R6_2.5.1 DelayedArray_0.14.1
#> [40] bitops_1.0-6 spatstat.utils_2.2-0 assertthat_0.2.1
#> [43] promises_1.1.1 scales_1.2.0 gtable_0.3.0
#> [46] globals_0.13.0 goftest_1.2-2 rlang_1.0.3
#> [49] RcppRoll_0.3.0 splines_4.0.2 rstatix_0.6.0
#> [52] lazyeval_0.2.2 spatstat.geom_2.3-0 broom_0.8.0
#> [55] yaml_2.3.5 reshape2_1.4.4 abind_1.4-5
#> [58] backports_1.4.1 httpuv_1.5.4 tools_4.0.2
#> [61] ellipsis_0.3.2 spatstat.core_2.3-1 jquerylib_0.1.4
#> [64] RColorBrewer_1.1-3 ggridges_0.5.2 Rcpp_1.0.8
#> [67] plyr_1.8.6 progress_1.2.2 zlibbioc_1.34.0
#> [70] purrr_0.3.4 RCurl_1.98-1.2 prettyunits_1.1.1
#> [73] rpart_4.1-15 openssl_2.0.2 deldir_1.0-6
#> [ reached getOption("max.print") -- omitted 88 entries ]